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ABSTRA CT 

A constrained iterative image restoration method is applied to 
multichannel diffraction-limited imagery. This method is based on the 
Gerchberg-Papoulls algorithm utilizing incomplete information and partial 
constraints. The procedure is described using the orthogonal projection 
operators which project onto two prescribed subspaces iteratively. Sprue of 
its properties and limitations are also presented. The selection of 
apropriate constraints was emphasized in a practical application. 
Multichannel microwave images, each having different spatial resolution, were 
restored to a common highest resolution to demonstrate the effectiveness of 
the method. Both noise-free and noisy images were used in this investigation. 
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I. Background 

Multichannel microwave radiometers on the Seasat and Nimbus 7 satellites 
offer a quantitative method for measuring geophysical parameters over the 
ocean. The emissivity of the ocean surface is low and varies predictably with 
wind speed; it thus provides a good background for observing precipitation. 
The theory and initial validation of this concept was given by Wilheit, et al. 

Cl]. 

Recently Olson [2] employed a radiative transfer model to simulate the 
polarized brightness temperatures that a Scanning Multichannel Microwave 
Radiometer (SMMR) would measure from hurricanes over sea surfaces at several 
frequencies (6.6, 10.7, 18.0, 21.0, and 37.0 GHz each with two 
polarizations). These brightness temperatures depend upon the rainfall rates, 
rain column height, and the emissivity of the wind roughened sea surface. It 
became evident that the 37 GHz channel is most sensitive to the height of the 
freezing level, whereas the 10.7 or 18 GHz channels are sensitive to changes 
in rainfall rates when those rates are less than 10 to 20 mm/hr. The 6.6 GHz 
channel provides rain information when the rainfall rate exceeds these 
values. The information content of each channel is a variable function of 
rainfall rate, rain column height, and emissivity of the sea surface, and the 
dependencies are generally nonlinear. A piecewise-1 inear regression algorithm 
has been applied to the synthetic data in the manner discussed by Smith and 
Woolf [3] to infer rainfall rates. The regression method employs data from 
eight of the SMMR channels. 

Unfortunately, the size of the antenna of the SMMR on Nimbus-7 imposes a 
diffraction limit on the sensor's angular resolution such that the relative 
angular response, d, of the radiometer is a function of the channel 
frequency. The antenna response function can be approximated by the 


t diffraction pattern of a circularly symmetric aperture with uniform 
illumination. The normalized radiation pattern, or the response function d, 
is of the form 


r'/O . (fea si n0)-i 2 

d < e > * IrVsTnd 


( 1 ) 


where a is the antenna radius, o is the angular deviation from the antenna 
centerline, k * 2ir/\ is the wave number and ( • ) is the first-order Bessel 
function. SMMR channels at different frequencies therefore have different 
footprint sizes, where the term footprint, or instantaneous field of view, is 
the most frequently used definition of the sputial resolution of a satellite 
radiometer. The foot prints corresponding to the half-power beamwidths of 
each channel are shown in Table 1. It is difficult to apply the regression 
algorithm unambiguously to real SMMR data, because each channel measures 
radiation from a footprint which may contain differing amounts of rain. This 
study overcomes the diffraction limitation imposed on spatial resolution by 
means of a constrained iterative restoration algorithm. 
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PARAMETER 

CHANNEL 

1 

2 

3 

4 

5 

WAVELENGTH (CM) 

4.54 

2-8 

1.66 

1-36 

0-81 

FREQUENCY (GHz) 

6.6 

10-7 

18-0 

21-0 

37.0 

DYNAMIC RANGE (° K) 

10 - 330 

FOOTPRINT SIZE (Km 2 ) 
(PICTURE ELEMENT) 

148*95 

91*59 

55*41 

44*30 

27*18 


Table 1. 
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II, Iterative Image Restoration - An Introduction 

The problem of image degradation can be stated mathematically as 
follows: the observed image, b(x,y), is expressed as 

b(x,y) = *(x,y) * d(x,y) + n(x,y), (2) 


where 

d(x,y) is the point spread response function of the degradation, 

*(x,y) is the ideal image, 
n(x,y) is random noise, and 
* denotes a convolution operator 

i 

In the frequency domain, we can also represent the spectrum of the 
degraded image by 


B(w x ,w y ) = L(u> x ,w y ) • D(w x ,w y ) + N(o> x ,u> y ) (3) 

where B, L, D, N are the Fourier transformations of b, a, d, and n, 
respectively. 

One approach to the image restoration problem is based on the concept of 
inverse filtering, in which the degradation point spread function is inverted 
to obtain a restored image. This approach is by solving for either an 
estimate of the ideal image z or its transform L. 

There have been a number of methods developed in recent years to restore 
degraded images, such as Wiener filtering, parametric estimation filtering, 
and pseudoinverse spatial restoration; see Andrews and Hunt [4], Wiener 
filtering is perhaps the most popular linear restoration method in which the 
mean square error between the ideal image and the restored image is 
minimized. Other commonly used restoration methods include the power spectrum 
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equalization method that sets the power spectrum of a blurred image equal to 
that of the original image, and the minimum a posteriori density method that 
uses Bayes theorem to optimize the restoration. 

The major task in linear restoration is to invert the degradation 
equation which is modeled either by a superposition integral or by a vector- 
space matrix equation. The existence of a unique inverse of the degradation 
point spread function, d(x,y), is required to derive £(x,y) in the presence of 
noise, n(x,y). If the inverse transformation of d(x,y) does not exist, then 
the problem is said to be singular. In those cases for which the inverse 
exists but is not unique, multiple solutions result. However, even if a 
unique inverse of d exists, it may be ill-conditioned. It can be shown that 
an arbitrary small perturbation in b can result in a large perturbation 'in the 
solution, i . Namely, the noise n(x,y) will lead to an undesired solution if 
the inverse of d is ill-conditioned. 

For these reasons, an alternative approach that is based upon iterative 
processing (successive approximations) to image restoration is of interest. 
This approach has the advantage that it is not necessary to determine the 
inverse of the point spread function. It has been demonstrated in a number of 
studies and applications that the sequence of approximations converges to a 
unique solution if constraints based upon known properties of the desired 
solution are incorporated into the iterative procedure. The use of iterative 
methods provides the flexibility of mixing constraints and distortions which 
is essential in solving many practical restoration problems. It is also 
evident from past work that this approach may provide a relatively efficient 
solution for linear spatially varying and nonlinear degradation problems (see 
Schafer, et al. [8] for a survey of iterative restoration algorithms). 
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For diffraction-limited imaging, the restoration of the image beyond the 
diffraction limit r annot be accomplished by filters that are based upon an 
inversion of convolution equation. The removal of noise and distortion is 
obtained at the expense of resolution. However, the use of iterative 
restoration techniques may provide superresolution in the presence of noise, 
i Such techniques make use of an incomplete a priori knowledge that we may have 
about the original image and the imaging system in order to extrapolate the 
bandlimited image spectrum beyond its diffraction limit. The incorporation of 
pnor knowledge of the image signal in the restoration process may be 
considered as a constraint operator which operates on the image signal 
iteratively to correct for imaging distortions. 

An iterative restoration technique using a spectral extrapolation 
algorithm was developed by Gerchberg [5] for the reconstruction of the signal 
of an object with a known finite duration or extent. In this method the a_ 
priori constraint operates on the signal's magnitude and the magnitude of its 
Fourier transform. Papoulis [6] proposed the same algorithm assuming a 
bandlimited signal and emphasizing signal extrapolation. Both of these 
methods are based on the fact that a finite object has an analytic spectrum. 
Analyticity implies that knowledge of only part of the spectrum is sufficient 
to uniquely determine the remainder of the spectrum. Hence, the complete 
spectrum may be derived from a diffraction-limited image although the 
information content may be incomplete. 

The Gerchberg-Papoul is algorithm has been studied and applied by a number 
of researchers. Sabri and Steenaart [7] derived a fast bandlimited signal 
extrapolation method using a matrix formulation. Another similar approach to 
the extrapolation of a bandlimited signal was proposed by Cadzow [16]. 
Rushforth and Frost [9] introduced a modified version of the Gerchberg- 
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Papoulis algorithm that provides additional flexibility and reduces the 
sensitivity to noise, Youla [10] generalized the algorithm using an 
alternating orthogonal projection method in Hilbert space and developed 
general conditions for convergence. Stack, Cahana, and Webb [11] formulated 
the restoration of finite-energy objects from two projections, Cahana and 
Stark [12] extended the algorithm to produce a faster convergence rate. 
Howard [13] formulated the algorithm as the solution of a system of linear 
equations containing Fourier components. Maeda and Murata [14] attempted to 
suppress noise in the restoration algorithm by preprocessing the image by a 
conventional spatial frequency filter. Fienup [15] proposed a version of the 
iterative algorithm for reconstruction using the magnitude of the Fourier 

i 

transform of a signal with positivity constraints. More recently, Hayes, Lim, 

and Oppenheim [17], [18] have investigated iterative procedures for the 

reconstruction of signals with finite extent from either the phase or the 

magnitude of the Fourier transform. All of these methods fit into the general 

framework of a general iterative restoration algorithm, with the Gerchberg- 

Papoulis algorithm being perhaps the most widely know algorithm of this type. 

In the following section, we will discuss the mathematical framework, 

convergence criteria, and the effect of noise in the general Gerchberg- 

« 

Papoulis algorithm. Some of the mathematical formalism and notation used by 
Youla [10] and Stark et al . [11] will be used here to describe the iterative 
restoration procedure. Experimental results obtained by applying a modified 
Garchberg-Papoul is algorithm to the restoration of multichannel microwave 
imagery will be presented. 
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_III> The Constrained Iterative Restoration >rithm 

A. The Resoration Algorithm 

Consider the Hilbert space H with elements f, g, h, x, etc., consisting 
of all L 2 square-integrable, Fourier-transformable functions, with inner 
product 

(f>9) s / f(x) • g(x) dx (4) 

-00 

and norm 

It fit B [ / [f(x)] 2 dx] 1(/ 2 . (5) 

-00 J 

Let and tf 2 be two subspaces of H. The angle between the two subspaces 
is defined by the expression 

1 

^(H- , H 0 ) = cos ( sup 

feH, 
geMg 

The angle ip=0 occurs wnen tfj and W 2 are 
each f e K we have a unique decomposition 

f = g + h, 

£ 

where g is the projection of f onto the subspace V and h is the projection of 
f onto the orthogonal complement of v> denoted by x ?. Since g and h are 
mutually orthogonal, the inner product (g,h) = 0. 

Two linear operators P(») and Q(») are defined by the rules 

g = P(f) and h = Q(f), (8) 

* 

respectively. p and Q are orthogonal operators which project the 
function f e H onto the orthogonal subspaces V and x 9, respectively. 


(Mil \ 

tl - nail / ' 


fll ' II g 


( 6 ) 


orthogonal to each other. For 


( 7 ) 
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Let us associate these two projection operators with the transformation 
of image signals between the spatial domain and the frequency domain via the 
Fourier transformation pair. Let p (•) and Q (•) denote the projection 
operators in the spatial domain which project onto p and j. p , 
respectively. Similarly, let ( • ) and Q b (‘) be the operators in the 
frequency domain which project onto p b and i <p b , respectively. 

For example, if a portion of the ideal image f(x) is known a priori (say, 
we have an unknown object in a background of known intensity), the unknown 
portion of f(x) can be expressed using the spatial projection operator, P a , as 

M(x) • f(x) « P a (f) , (9) 

while the known portion of the image f(x) is expressed using the projection 
operator Q a as t 

(1 - M(x) ) • f (x ) = Q fl (f) , (10) 

where 

[ 1, xe unknown extent 

0, x e known extent , 

This pair of spatial projection operators decomposes the image signal f 
into two parts: the unknown portion of the signal P a (f), and the known portion 
Q a (f), where f = P a (f) + Q a (f) and (P a (f), Q a (f)) = 0. 

By dividing the image spectrum into two spectral bands at u> c , that is, a 
low-pass spectrum with Jwj < u> c and a high-pass spectrum with jw| _> w , the 
image signal f is decomposed by two projections f = P b (f) + Q b (f) with 
(Pb(0> Qb ( ^ ) ) = °* operators P b ( # ) and Q b ( # ) are defined by the 

following equations: 
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«w<aaa«nMr ■ w wi— iw.'** 1 •'■*"• 


*F" 1 CH(w) • F(ui)3 ■= P b (f) (11) 

T _1 C(1 - M(«)) • F(»)3 - O b (f), (12) 

is the Fourier transform of f(x)> 
dene! •( J inverse Fourier transformation, 
are projection operators onto the low-pass and high-pass 
image subspaces, respectively, and 

Hw-r 1, for w< “ c 

■* 0, for |w| J> w . 

i 

The iterative image restoration problem can be stated as follows: given 

the projections of the image signal on prescribed subspaces, can one derive 
the original image f (x ) ? It has been proved by Stark, et al. [11] that the 
restoration is possible under certain conditions if partial information from 
the subspaces is given. 

To be more precise, let us consider an example: If we are given two 

projections of f(x), say g - Q a (f) and h = P b (f) as defined in (10) and (11), 
respectively, then a restoration is possible. In this case, g - Q a (f) is 
considered to be the portion of the image data which vanishes outside some 
known extent, and h = P^f) is considered to be the low-pass signal of a 
bandlimited image. The low-pass spectrum, f[Pj ) (f)], is the truncated spectrum 
of the image f. This low-pass spectrum can be obtained from the fact that the 
low-frequency components of the degraded image have not been changed, or the 
knowledge of the degradation function is known. This example is equivalent to 
the extrapolation problem described by Gerchberg [5] and by Papoulis [6]. 


where 

F(«) 

P b and Q b 


Before starting the iterative restoration algorithm, wo first define the 

00 


notation of iterative projection, (P Q b ) ' (f ), explicitly by 
(P a Q b ) (k, (0 ■ P a <V WV”W f >)))...). 


(13) 


where both P a and (V oper '■ times. For example, 


(P a *. (0 - f 

(P a C. '(f) » P a (O b (f)) (14) 

( p a f ;,! !2) ( f ) - P a (Q b (P a (Q b (f)))) • 


Then, the iterative method is stated as follows: 

i 

Step l: Obtain the signal P b ( f ) and Q a (f), based on available a pri.‘ri 

information, 

Step 2: Let k=l, and let the first guess of the iteration process be 

f(1) ■ P a ( p b ( f )) + Q a (f). (15) 

Note that fW is the enhanced image after k-1 iterations. 

Step 3: Correct the error by imposing constraints in the spectral domain. 

f‘ ■ f( k) - P b (f ( . k> ) + P b (f) (16) 

Step 4: Correct the error by imposing constraints in the spatial domain, 

f( k+1 > » f> . Q a (f * ) + Q a (f) (17) 

Step 5: Let k = k + 1. If a satisfying result is obtained, stop. 

Otherwise, go to step 3. 

i 

Figure 1 shows the block diagram of the restoration method. 
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a priori Spatial 
Knowledge Q a (f) 



In Hoi Guess 


f < 0 »P a (P b (f))+Q a (f) 










t 

The procedure begins with an initial guess as in (15) by performing the 
operation P fi on the signal P b (f) to remove all the signal from regions where £ 
priori data are available. Adding Q a (f) then restores the correct values back 
into those regions. Both the correct va’ues and their corresponding spatial 
locations are assumed to be known. This initial guess, f(’), hi cne result of 
a one-step Gerchberg-Papoulis algorithm, which performs error correction on 
the low-pass filtered image, P b (f), using incomplete information. 

In step 3, spectral constraints are imposed according to (16) where 
incorrect low frequency components at the k th iteration, P b (f( k )), are 
replaced with spectral information P b (f) that is assumed to be known. In step 
4, spatial constraints are imposed according to (17); the resulting signal f 1 

i 

from (16) is corrected by removing Q a ( f ' ) and replaced by the known spatial 
information Q a (f), 

It is easy to show that Equations (16) and (17) are functions of f^" 1 ), 
the previous iteration result, and the initial guess as in Eq. (15). 

Substituting Eq. (16) into Eq. (17), and letting k = 2,3,..,, we obtain 

+ f<1> 

f (3) = P a (Q b (f (2) )) + f (l) (18) 

• * 

• 

• 4 • 

= p a (Q b (^ ( k ~ l h) + 

Solving these equations recursively, a simple equation is obtained for 
the k^ image iterate: 
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( k ) 


k-1 


(vi) 


f'" 1 - ); (p a Q b ) uy ( f - (p a Q b )(f» 

j-Q 

- f - (P a Q b ) (k) (f) 


(19! 


The reconstructed image f( k ) can thus be considered as the difference 
between the ideal image f and the k~times projected version of the ideal 
image, (P g Q b )^)(f). in other words, f(^) depends on the ideal image function 
f, the projection operation P a Q b , and the number of iterations. If f( k ) is an 
estimate of f, then e( k ) s ( p a%)^^ f ) is the error in the kth irTia 9 e derate. 
Since Q b and P a are projection operators, the norm (or power) of the 

projection onto the subspace is less than or equal to that of the mapping 

signal, that is, 

» 

1 

llfll 2 HQ b (f)n, 

lif ll _> BP (f )ll , and (20) 


Hf« > ll(P a Q b )(f)H • 

The error term in the (k+l) st iterate, e^ + l), can be expressed in terms of 
e^ k ), the error in the k^ iterate by 

e (k+1) - (P s Q b ) (k+1) (f) (21) 

■ p a%(< p aQ b f (k) < f >> 

= (P a Q b )(e (k) ). 

with 


1 e ( k + 1 } H < n(P a Q b )(e' k ^)i 


(22) 
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This implies that the error is strictly monotonical ly decreasing, and it has 
been proved [10] that recursion Eq. (19) converges to f as k approaches 
infinity., if and only if v n 1 1> b « { } and y[v a , i 7> b ) > 0. 

In practical situations the iterative process converges to some 
constant not equal to f for some finite k, in which (P a Q b )^(f) = 

(P a Qb ) (f)< This is mainly due to noise generated by imaging, 
quantization error, an inaccurate estimate of the degradation function, etc. 
One important noise source that will further affect the performance of the 
restoration is the insufficiency of a priori constraints. In this situation 
the prescribed subspaces do not provide sufficient information for the 
restoration when \'j(’P a , 1 *P b ) is very small. Such a situation can be improved 

i 

by incorporating more constraints into the procedure, 

* 

B. Restoration with Noisy Constraints 

Insufficient or inaccurate constraints imposed on the restoration 
procedure will lead to undesirable results even if the iteration process 
converges. For example, if we have some uncertainties in the projection 
operators which project onto subspaces v and jl V, the extrapolated signal will 
inevitably contain some error, and this error will amplify in the iteration 
process. In most cases, noisy constraints 'are the result of wrong assumptions 
in the physical model of the problem, or from noisy image signals used to 
constrain the data. 

We now demonstrate the effect of noisy constraints on image resolution 
restoration. We use two projections of f(x), say g' and h', where 
9 1 “ (Q a (0 + n s ) ar| d h* = (P b (f) + n f ), to represent the amount of additive 
noise embedded in the partial information from the two prescribed subspaces. 
In the case of bandlimited image extrapolation, n s and n^ represent erroneous 
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constraints in the spatial and frequency domains, respectivly. The initial 
guess of the iterative procedure becomes 

f (1) - P a (P b (f) + n f ) + (Q a (f) + n s ). (23) 

By substituting f = P fl (f) + Q.(f) and f « P b (f) + ( f ) into (23), we have 

f(1) ' p a ( p b ( f )) + f - p a ( f ) + p a< n f) + n s t 24 ) 

s [f “ P a( Q b{ f )^ + [P a^ n f^ + n s^* 

Rewriting (16) and (17) using the noisy constraints, we obtained the recursion 
formula 

f (k) ■ Y ( p a Q b ) (J) (f - P a (Q b (f)) + P a Cn f ) + n s ) (1-5) 

■ f - ( p a0b) (k) ( f ) + ^ ( p A> (J) ( P a<"f> + n s ) 

=f-e< k >. 

The error e^) i n the restored image f^) consists of two terms. The 

first term (P-Qu)^(f) decreases monotonical ly while the second term 
k-i m . 

I (P a Qb ) '( p a ( n f ) + n g ) increases monotonical ly as the number of iterations 
k increases. One must therefore consider the tradeoff between these two error 
terms in deciding how many iterations should be applied. If the image signal 
is much stronger than the noise signal, the iterative procedure will converge 
within a reasonable number of iterations. For the case of extremely noisy 
constraints, the second error term becomes dominant. Thus, fewer iterations 
will provide satisfying results while more iterations will lead to divergence. 

In one experiment a square object in a square background field, both with 
constant signal levels, was used to demonstrate the tradeoff between 
performance and the number of iterations. The test image was first degraded 
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by a low-pass filter. The bandllmlted image was then restored by the above- 
described constrained iterative algorithm. The cutoff frequency of the image 
signal w«s assumed to be known and it was used as a constraint in the 
frequency domain. The physical extent of the square object and the signal 
from the background were considered as partial a priori spatial information. 
To generate the effects of noisy constraints, independent noise was added to 
the spatial constraint. This was done by corrupting the signal of the 
background with additive zero-mean white Gaussian noise. Noise levels with 
different standard deviations (0, 5, 10, 20, 30, and 40) were used. The 
restored images were evaluated by measuring the mean-square error between the 
ideal image and the restored image. Figure 2 shows the resultant curves 

t 

representing the performance of the algorithm at different noise levels. At 
low noise levels, the procedure converges in a reasonable number of iterations 
to a low error rate. This implies that the bandlimited image has been 
successfully restored. With higher noise levels in the constrained data, the 
second term of in (25) becomes dominant and the algorithm performance 

deteriorates as k increases. This experiment confirms that the restoration 


results are sensitive to the number of iterations when noisy constraints are 


used. 



Figure 2. Performance of the algorithm with noisy constraints 
at different noise levels. 
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IV. Restoration of Multichannel Microwave Images 

The image restoration procedure described above has been aplied to 
overcome the diffraction-limited imagery obtained from the SMMR, The goal is 
to restore channels of low spatial resolution to a common high resolution. 
The low frequency bandlimited images are obtained from the 21.0, 18.0, 10.7, 
and 6.6 GHz channels, and the common high resolution is that of the 37 GHz 
channel s. 

The procedure described in Section III was implemented with special 
features, which accommodates the multichannel SMMR imagery of hurricanes. The 
procedure is summarized in Figure 3. Details of the partial constraints will 
be discussed in the following sections. 

i 

A. Spatial Constraints, Q a (f) 

Spatial constraints are based on the estimated extent of rain areas, 
derived from the 37 GHz and infrared images, the upper and lower bounds of the 
measured brightness temperatures, and some physical attributes of hurricanes. 

We initially utilized polarization information at 37 GHz to provide an a_ 
priori estimate of the spatial extent of model hurricane raincells. While 
microwave emission from the ocean surface is highly polarized, the emission 
from raincells is nearly unpolarized. Therefore, a polarization threshold was 
applied to the 37 GHz antenna temperatures to identify areas of significant 
rainfall (rainfall rate > 2 mm/hour). 

A simple physical model was then employed to relate the 37 GHz antenna 
temperatures outside the rain areas to antenna temperatures at the lower 
microwave frequencies. In these regions, the microwave antenna temperature, 
T b , at frequency v and in polarization p can be expressed as 
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Figure 3. The procedure to restore low frequency SMMR images 

(6.6, 10.7, and 18.0 GHz) to a high common resolution 
(that of the 37 GHz channels). 











V v > P) ■ (1 - r(v, p)) • T ir «(,) 


( 21 ) 


+ T a (v) • (1 - n(v)) • (1 + r(v,p) • n(v) ) 

+ T bb • • r ( v >p)> 

where T A (v) - Tj(v) + T 2 (v) • (1 - n(v)) 

is the effective emitting temperature of the atmosphere. = 2.7 °K is the 
cosmic microwave background, and T IR is the infrared ocean surface 
temperature, which can be obtained operational ly from geostationary 
satellites. The constants T^(v) and Tg(v) are obtained empirically. The 

atmospheric transmittance is JI(v), and the ocean surface reflectivity, 
r(v,p), can be expressed [19] 1 

V 

r(v, p) = a(v, p) - 3(v) • U 2Q , (22) 

where U 20 is the wind speed at 20 m height, and a(v, p) and e ( v) are empirical 
constants. The three terms on the right-hand side of (21) represent surface 
emission, atmospheric emission, and the small contribution from the cosmic 
background, respectively. 

From (21) and (22) it may be noted that each antenna temperature depends 
only upon the two unknown quantities n(v)„and U 2Q . Thus, antenna temperature 
measurements in the two polarizations at 37 GHz provide enough information to 
solve simultaneoul sy for n(37 GHz) and U 2Q . Since the transmittance at 37 GHz 
is related empirically to transmittances at the lower frequency channels, 
Tg (21.0, p), Tg (18.0, p), Tg (10.7, p) and Tg (6.6, p) can also be computed. 

In areas where the rainfall rate did not exceed ~ 2mm/hour the above 
method (Eqs.-(21) and (22)) was utilized to recreate imagery in the low 
frequency channels with 37 GHz resolution. This "known" portion of each low 
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frequency image was then used as a partial constraint, or incomplete 
information, to extend the resolution in the heavier rain areas, 

B. Spectral Constraints, Ph(f) 

The spectral constraints used in the restoration algorithm are based on a 
knowledge of the highest cutoff frequency (defined by the 37 6Hz resolution), 
and a knowledge of the degradation point spread function. The SMMR antenna 
aperture illumination function is considered to be the degradation point 
spread function of the imaging process. 

In this initial study, a Bessel -type point spread function v/as used to 
describe the passive microwave radiometer response; see Eq. (1). This 
response function is a simplified model of the actual SMMR antenna aperture 
illumination that allows us to assume a flat response function at 'the low 
frequency portion of the spectrum. During the iteration, the DC and a few low 
frequency components (both the phase and amplitude) of the image spectrum were 
kept intact, while frequency components above 37 GHz were removed. 

This simple response function (1), with the assumption of a uniformly 
illuminated aperture, does not agree in general with the actual SMMR antenna 
aperture illumination function for which one must assume non-uniform 
distribution. In order to apply the restoration algorithm to real microwave 
image data, a more realistic degradation point spread function will be needed. 

One can empirically estimate the amplitude and phase distribution 
functions of the SMMR instrument. One can also model the antenna response 
through analytical methods. The most commonly used phase distributions across 
the aperture are linear and quadratic functions, where linear-phase 
illumination is the basic principle behind an electronic scanning antenna, and 
quadratic-phase illumination functions are used primarily to effect far-field 
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conditions in the Fresnal region of the diffracted field. The amplitude 
illumination is governed by antenna parameters such as the antenna 
directivity, the side-lobe energy level, and the half-power beamwi dth . The 
amplitude illumination function of a particular antenna can be modeled by 
matching the antenna parameters with those produced by selected distribution 
functions, such as the cosine functions, the parabolic functions, etc. 

Given an estimate of the point spread function -- both in phase and 
amplitude — we will be able to constrain the spectrum during the iterative 
restoration process and apply the procedure to restore the 6.6, 10.7, 18.0, 
and 21.0 GHz SMMR channels to a common, high resolution (that of the 37 GHz 
channels). 

, i 

, * 

C. Experiments and Results 

A set of synthetic hurricane images, each consisting of antenna 
temperatures in a 16 x 16 pixel image field, were created. Noise-free antenna 
temperatures at 37.0, 18.0, 10.7, and 6.6 GHz were generated using the 
approximate radiative transfer model of Olson [3], assuming an ocean surface 
background (see Figure 4). After the restoration, the low frequency channels 
were enhanced to a resolution compatible with the 37 GHz channel. The 
enhanced images are shown in Figure 5. . 

Figure 6a shows the cross-sections of rain cells in the 6.6 GHz synthetic 
hurricane image before and after restoration. The degraded 6.6 GHz image 

resolution has been enhanced to a large degree after a few iterations. A 
synthetic image of the entire model hurricane that would be measured at 37 
GHz, the original degraded 6.6 GHz image, and the enhanced 6.6 GHz model 
hurricane image are shown in Figs. 6b, 6c and 6d respectively. 
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Figure 4. The observed 6.6, 

10.7, and 18.0 GHz 
images at the parallel 
polarization. 
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Figure 6. Enhancement of the resolution of a synthetic 6.6 GHz image 

of a hurricane. A detailed cross-section plot is shown in (a). 
This enhancement algorithm uses the known extent of the image 
shown in (b) to enhance the degraded 6.6 GHz image shown in (d). 
The derived image at 37 GHz resolution is shown in (c). 
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Data distorted by noise has rendered It difficult to continue the 
spectrum beyond the original diffraction limitation. The restoration 
procedure was applied to a set of noisy Images to determine Its ability to 
handle noisy data. In one example, additive white noise with an rms value of 
4° K was added to both the 37 GHz and 6.6 GHz synthetic linages. Within a few 
iterations, we obtained the restored image scan as shown in Fig. /, The 6.6 
GHz noisy image Is almost completely restored to the optimal resolution. 

V. Conclusions 

An iterative resolution restoration method for restoring multichannel 
diffraction-limited imagery has been described. It is based upon the 
Gerchberg-Papoul 1 s algorithm using incomplete information and partial 
constraints in both the image space and the Fourier space. The procedure was 
presented using the orthogonal projection formulation proposed by Voula [10] 
where projection operators project onto two prescribed subspaces. The 
projection operators are defined by incorporating a priori information of the 
imaging system and the nature of the observed object. The tradeoff between 
error and convergence rates has also been investigated. 

One of the focuses of this research was on the multichannel microwave 

t 

image restoration problem, with special attention paid to the selection of 
appropriate constraints for the iteration procedure. It was demonstrated that 
the constraints control the performance and rate of convergence of the 
restoration. Using the procedure, a set of synthetic hurricane images was 
restored to a common resolution. 

This application has shown the effectiveness of this iterative procedure 
in restoring the spatial resolution of low frequency channels. An obvious 
next step is to apply this procedure to real multichannel data. 
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